#ifndef _STABLEINSSOLVER2D_H_
#define _STABLEINSSOLVER2D_H_

#include <cmath>
#include <algorithm>


class StableINSSolver2D {

private:

  int domainSizeX, domainSizeY;

  float hx, hy;

  int gridSizeX, gridSizeY, gridSize;

  float *velocityX, *velocityY;

  float *tmpVelocityX, *tmpVelocityY;

  bool *boundary;

  float *tmpGrid;

public:

  enum BoundaryCondition { BC_DIRICHLET, BC_DIRICHLET_ZERO, BC_NEUMANN_ZERO };


  StableINSSolver2D();

  ~StableINSSolver2D();

  const float *getVelocityX();

  const float *getVelocityY();

  void resetSystem(
      int domainSizeX, int domainSizeY, float gridSpacing,
      const float *velocityX, const float *velocityY,
      const bool *boundary);

  void updateVelocity(
      float *sourceX, float *sourceY,
      float viscosity, float vorticity, float dt,
      BoundaryCondition bc = BC_DIRICHLET_ZERO);

  void transportScalarField(
      float *field, float *source,
      float diffusion, float dissipation, float dt,
      BoundaryCondition bc = BC_NEUMANN_ZERO);

private:

  void copyBoundaryValues(float *target, const float *source);

  void addSource(
      float *target,
      const float *field, const float *source,
      float dt);

  void addVorticity(
      float *u, float *v,
      float vorticity, float dt,
      BoundaryCondition bc);

  void diffuse(
      float *field, const float *field0,
      float factor, float dt,
      BoundaryCondition bc);

  void advect(float *field, const float *field0, float dt, BoundaryCondition bc);

  void projectVelocity(
      float *pressure, float *divergence, float dt,
      BoundaryCondition bc);

  void dissipate(float *field, float dissipation, float dt);

  void solveLinearSystem(
      float *x, const float *x0,
      float a, float ai, float aj,
      BoundaryCondition bc);

  int idx(int i, int j);

  float value(
      const float *x,
      int i0, int j0,
      int i1, int j1,
      BoundaryCondition bc);

  void freeData();
};

#endif /* _STABLEINSSOLVER2D_H_ */

